!
! astro transforms
!
!
! Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

module astrotrafo

  implicit none
  
  integer, parameter, private :: dbl = selected_real_kind(15)
  real(dbl), parameter, private :: rad = 57.295779513082322865_dbl

  ! Transformation and projection:
  !
  ! type ... type of projection (GNOMONIC,..)
  !      ... mostly affine transformation ...
  ! scale ... scale [deg/pixel]
  ! pa    ... position angle [deg]
  ! acen ... center projection in alpha [deg]
  ! dcen ... center projection in delta [deg]
  ! xcen ... horizontal center projection [pixel]
  ! ycen ... vertical center projection [pixel]

  type TypeProjection
     character(len=80) :: type   
     real(dbl) :: acen,dcen,xcen,ycen,scale,pa
  end type TypeProjection
  
  private :: ptrafo,invptrafo

contains

  subroutine trafo(t, a, d, u,v)

    use projections

    type(TypeProjection), intent(in) :: t
    real(dbl), dimension(:), intent(in) :: a,d
    real(dbl), dimension(:), intent(out) :: u,v

    if( t%type == "GNOMONIC" ) then
       call ptrafo(gnomond,t,a,d,u,v)
    else if( t%type == "IDENTITY" ) then
       call ptrafo(identity,t,a,d,u,v)
    else
       write(*,*) "Required projection not implemented yet:",t%type
       stop 666
    end if

  end subroutine trafo

  subroutine invtrafo(t,u,v,a,d)

    use projections

    type(TypeProjection), intent(in) :: t
    real(dbl), dimension(:), intent(in) :: u,v
    real(dbl), dimension(:), intent(out) :: a,d

    if( t%type == "GNOMONIC" ) then
       call invptrafo(invgnomond,t,u,v,a,d)
    else if( t%type == "IDENTITY" ) then
       call invptrafo(invidentity,t,u,v,a,d)
    else
       write(*,*) "Required projection not implemented yet:",t%type
       stop 666
    end if

  end subroutine invtrafo

  subroutine ptrafo(proj, t, a, d, u,v)

    interface
       subroutine proj(a,d,a0,d0,x,y)
         implicit none
         integer, parameter :: dp = selected_real_kind(15)
         real(dp),intent(in) :: a,d,a0,d0
         real(dp),intent(out) :: x,y
       end subroutine proj
    end interface
    
    type(TypeProjection), intent(in) :: t
    real(dbl), dimension(:),intent(in) :: a,d
    real(dbl), dimension(:),intent(out) :: u,v
    real(dbl) :: x,y,aa,bb
    integer :: i

    if( size(a) /= size(d) .or. size(d) /= size(u) .or. size(u) /= size(v) )then
       write(*,*) 'Astro Trafo: dimension check failed.'
       stop 31
    end if

    aa = t%scale*cos(t%pa/rad)
    bb = t%scale*sin(t%pa/rad)

    do i = 1, size(a)
       call proj(a(i),d(i),t%acen,t%dcen,x,y)
       u(i) = t%xcen + aa*x - bb*y
       v(i) = t%ycen + bb*x + aa*y
    end do

  end subroutine ptrafo


  subroutine invptrafo(invproj, t, u, v, a,d)

    interface
       subroutine invproj(a,d,a0,d0,x,y)
         implicit none
         integer, parameter :: dp = selected_real_kind(15)
         real(dp),intent(in) :: a,d,a0,d0
         real(dp),intent(out) :: x,y
       end subroutine invproj
    end interface
    
    type(TypeProjection), intent(in) :: t
    real(dbl), dimension(:),intent(in) :: u,v
    real(dbl), dimension(:),intent(out) :: a,d
    real(dbl) :: x,y,xx,yy,aa,bb
    integer :: i

    if( size(a) /= size(d) .or. size(d) /= size(u) .or. size(u) /= size(v) )then
       write(*,*) 'Trafo: dimension check failed.'
       stop 31
    end if

    aa = cos(t%pa/rad)/t%scale
    bb = sin(t%pa/rad)/t%scale

    do i = 1, size(a)
       x = u(i) - t%xcen
       y = v(i) - t%ycen
       xx =   aa*x + bb*y
       yy = - bb*x + aa*y
       call invproj(xx,yy,t%acen,t%dcen,a(i),d(i))
    end do

  end subroutine invptrafo

end module astrotrafo
